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ABSTRACT 

There has been disagreement about whether cooling in protoplanetary disks 
can be sufficiently fast to induce the formation of gas giant protoplanets via grav- 
itational instabilities. Simulations by our own group and others indicate that this 
method of planet formation does not work for disks around young, low-mass stars 
inside several tens of AU, while simulations by other groups show fragmentation 
into protoplanetary clumps in this region. To allow direct comparison in hopes 
of isolating the cause of the differences, we here present a high resolution three- 
dimensional hydrodynamics simulation of a protoplanetary disk, where the disk 
model, initial perturbation, and simulation co nditio n s are essentially identical 



Bossl (120071 . hereafter B07). As 



to those used in a recent set of simulations by 
in earlier papers by the same author, B07 purports to show that cooling is fast 
enough to produce protoplanetary clumps. Here, we evolve the same B07 disk 
using an improved version of one of our own radiative schemes and find that 
the disk does not fragment in our code but instead quickly settles into a state 
with only low amplitude nonaxisymmetric structure, which persists for at least 
several outer disk rotations. We see no rapid radiative or convective cooling. We 
conclude that the differences in results are due to different treatments of regions 
at and above the disk photosphere, and we explain at least one way in which the 
scheme in B07 may lead to artificially fast cooling. 



Subject headings: accretion, accretion disks — hydrodynamics — instabilities - 
planets and satellites: formation — protoplanetary disks 
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Introduction 



The issue of Jupiter and Saturn's provenance, as well as that of the known extrasolar 
gas giant planets, remains a subject of intense study. Over the last decade, several 
computational groups, using an array of different techniques and disk models, have 



examined the disk instability mechanism, where gas-phase gravit ational instabiliti e s (G 
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20091 ). Researchers agree that GI's 



are triggered when disks are massive and cold and that, once GI's occur, cooling on time 
scales comparable t o the dynamic time is required for disk fragmentation (for a review, see 
urisen et al.ll2007l ). The questions are whether cooling by radiative transport with realistic 
opacities is rapid enough anywhere in protoplanetary disks to cause fragmentation and 
whether clumps that do for m remain physically intact, gravitationally bound entities before 



the gas disk dissipates (e.g., 



Haisch et al. 



20011 ) . There have been sharp disagreements on 



these points among the different groups attacking the problem, and so it is important to 
evaluate how and why different methodologies may lead to strikingly different conclusions. 

For moderate mass disks with a radial extent of tens of AU around solar- t ype st ars, 



simulations with radiative transport and realistic opacities presented in 



hereafter B07) support disk fragmentation, in partial agreement with 



but in severe d i sagre e ment with simulations by our own group (jBolev et al. 



Cai et al. 


2006 




2008 


Forgan et al. 


20091). 



Boley fe Durisen 



Boss! (12 007. 



Ma ver et al 



2006 



(2007) 



2007 



20081 ) and others (IStamatellos fe Whitworth 



2008 



2009). In B07, Boss considers a variety of factors that may account for 



the disagreement between his results and ours. We as well have explored some of the 
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issues described in B07, includi ng;: resolution (jPickett et al.ll2003t iBolev fc Durisenl 120081 ) 



radiative transport algorithms (Boley et al 



viscosity (jPickett &: Durisen 



20071 ). irradiation ( ICai et al. 



20081 ). and 



20071 ). In this paper, we focus on differences in the radiative 



schemes by using initial conditions and input physics that are as similar as possible to 
those described in B07 and by keeping as much of the numerical treatment as close to the 
techniques and conditions regularly used by Boss. By employing a different and, we think, 
better treatment of the radiation physics, we find, as we have before, that realistic cooling 
is not nearly strong enough to initiate clump formation. 

This paper is organized as follows. In Section 2, we describe our numerical methodology 
and initial axisymmetric equilibrium state. We describe the results for a simulation lasting 
five outer rotation periods in Section 3. In Section 4, we compare our results with those in 
B07 and try to isolate the causes for our differences. Section 5 is a brief summary. 



2. Computational Methodology 

2.1. 3-D Radiative Hydrodynamics 

We conducted our three-dimensional disk simulation using the CHYMERA 
(Computational HYdrodyn amics with Multip lE Radiation Algorithms) code, developed 



at Indiana University (e.g., 



Boley et al 



20071 ). CHYMERA uses an Eulerian scheme 



that is fully second-order in space and time to solve the equations of hydrodynamics and 
Poisson's equation on a cylindrical grid. The grid has a resolution of (256,5f2,64) in 
cylindrical coordinates (r, 0, z). The disk is initially located between radial cells 40 and 
202, corresponding to 4 and 20 AU, but is free to expand both radially and vertically. 
Hydrodynamic boundary conditions are outflow along the top, sides and inner hole of the 
grid; material from the disk that moves inside the inner radius at radial cell number 13 is 
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automatically added to the central star's mass. Equatorial plane symmetry is assumed. 



The pr imary difference be twee n the version of CHYM ERA used in this paper and the 



one used in 



Bolev et al. 


(2006 


) and 


Cai et al. 


(2006. 


2008) 



20081 ) is an improved treatm ent of the 



radia tive physics in the optically thin region (for more details and test results, see 
20091 ); Cai et al. 2010, in preparation). 



Zhu et al. 



Define r as the optical depth measured vertically downward from above, calculated for 
each computational cell using Rosseland mean opacities. Regions in which r > 2/3 are part 
of the disk interior, and we calcuate the radiant energy flow in all three directions using 
flux-limited diffusion with Rosseland mean opacities. On the other hand, we allow the 
atmospheric cells (in which r < 2/3) to cool radiatively via an optically thin LTE emissivity 
based on the local temperature but computed relative to radiative equilibrium. Thus, the 
volumetric cooling rate in any atmospheric cell is 

A = 4 P K(T)(aT 4 -it J), (1) 

where k(T) is the Planck mean opacity and J is the mean intensity, given by 

3d, 



nJ=^ h (r + 2/3) + aT ( 



4 

envi 



(2) 



which is a radiative equilibrium solution for an Eddington-like grey atmosphere including 
external irradiation with a downward flux at the upper boundary of uT^ nv . We couple the 
optically thick and thin regions using a boundary condition, which defines the r = 2/3 
boundary flux oT^ h at the photosphere for the interior diffusion approximation. 

This scheme effectively measures optically thin local heating and cooling relative to 
the radiative equilibrium represented by J and avoids the discont inuity at the photo sphere 



and the artificially reduced atmospheric temperatures reported in 



IT 



Bolev et all (I2006D for an 



earlier and cruder sche me. A similar r adiative scheme was successfully employed to study 



FU Orionis outbursts ( IZhu et al 



20091 ) 
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In B07, on the other hand, no radiative physics is done in the optically thin regions. 
Instead, the temperature above the photosphere is reset to that of the background 
temperature of an assumed thermal bath every time step, without computation of the 
detailed heating and cooling of optically thin cells. The boundary temperature for the 
radiative diffusion solution in the disk's optically thick interior is then taken to be the 
background temperature. 



he opacities and molecular weights in our simulations are described in iD'Alessio et al. 



( 120011 ) and are similar to those used in the Boss calculations. In the pres ent case, we ado pt 



a revised mean molecular weight of \x = 2.39 (lowered from the value in 



Cai et al 



20081 ) 



for the temperatur e and pr essure ranges in our disk simulation, close to the //-value used 



in B07 ( Boss 



200l|). As in 



Cai et al. 



(120081 ). we set the maximum particle radius in the 
dust opacity to one micron, compatible with B07. Note for clarity that we are not using 
the option in CHYME RA where the enti r e ^-direction is comput ed using radiative transfer 



along a single ray (e.g., 



Boley et al 



20071 : 



Boley &: Durisen 



2008|). 



Equation (2) includes irradiation as a downward black-body flux at r = from an 
envelope with temperature T env = 50 K. We do this to mimic the thermal bath approach 
in B07, where the optically thin regions are maintained at a fixed temperature, typically 
50 K. We point out that the effect of the thermal bath in B07 is to set a temperature 
boundary condition (BC) for the diffusive optically thick region, rather than a radiative 
flux BC to couple the thick and thin layers, as we do. As far as we know, the photospheric 
BC and the explicit radiative treatment of the r < 2/3 region in our code represent the 
only substantive differences between the physics included in the Boss simulations and our 
own. One minor diff erenc e in th e interior treatment is that we have enabled the flux limiter, 
which, according to iBossl (120011 ). should not change the outcome. On the numerical side, 



the B07 grid is spherical, not cylindrical, so the disk photosphere in B07 is located by a r 
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condition along spherical radii not the z-direction. 



2.2. Initial Model and Simulation Conditions 



The protoplanetary disk system is comprised of a star with one solar mass and a disk 
of mater ial in nearly Keplerian rotation with a mass equal to 0.091 M Q . The initial model 



uses the iBosd (11984] ) pressure equation of state (EOS) based on an initial temperature 
profile provided by Boss (2004, private communication), but we derive the specific internal 
energy density by using e = p/(j — 1), where 7 is the ratio of specific heats and e is the 
internal energy density. In the subsequent evolution, the EOS is assumed to be that of an 
ideal gas with 7 = 5/3. As in previous simulations, the disk begins in a marginally unstable 



state and is allowed to cool radiatively to insta bility. 



he in itial axisymmetric model is set 



Boss! (120031 ). but without infall. The latter 



up using the same analytic expressions given in 
is not a serious omission, because, in the B07 simulation, there is little mass in the inflow, 
and infall lasts only ~ free-fall time, about 0.2 ORP. We included this short accretion phase 



in preliminary simulations with the 



CaietaL 



(120081 ) version of CHYMERA and obtained 



essentially the same results reported here. We introduce the same initial perturbation used 



in B07 which preferentially seeds power at the per cent 



evel into ordered cos(mcb) structure 



1998 



20001 ) . We ran the simulation 



with m = 2 to 4 plus a small random component ( IBossI 
250,000 computation steps, equal to 5.0 ORPs. An ORP (outer rotation period) is the 
orbital period for material initially located at 20 AU, or about 90 yrs. The end of the 
simulation extends a little more than one ORP beyond the time shown in Figures 2 and 3 
in B07. Though material is allowed to expand off the grid (and accrete onto the central 
star), the disk itself loses less than 0.4 % of its original mass and angular momentum. 
The simulation was conducted on a dedicated HP Proliant DL385 G2 server at Lawrence 
University. 
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3. Results 

Figure 1 shows a three panel equatorial plane mass density contour map at three 
different times in the simulation. In each case, the density contours span about 7 orders 
of magnitude, with each contour level representing a factor of 2 change in density; filled 
contours indicate regions of the highest density (p > 10~ 10 gem -3 , following the convention 
of B07). The box for each image has dimensions of 40 AU x 40 AU. The top panel shows 
the model at the time nonaxisymmetric structure reaches its maximum strength, around t 
= 0.25 ORP. Note the strong low-order spiral disturbance (primarily four-armed) near the 
inner boundary, which is a result of the initial perturbation. In the middle panel, after 3.7 
ORP's, only a remnant of nonaxisymmetric structure is seen in the bulk of the disk. The 
highest density regions and largest amplitude disturbances are both located near the inner 
boundary of the grid, though some spiral features can be traced lightly in the outer regions. 
This is completely different from what is shown at the same evolutionary time for models 
H and TZ in Figures 2 and 3 of BO 7, namely, dense clumps and/or nascent clumping near 
10 AU. Over the final two ORP's of our simulation, as shown in the bottom panel of Figure 
1, our disk has evidently settled into a quasi-steady configuration, consisting primarily of a 
dense ring of material near the inner hole and a large but low-amplitude two-armed spiral, 
parts of which extend beyond the original 20 AU edge of the disk. Overall, the primary 
effect of the GIs is radial transport of mass. 



Unlike some of our previous simulations, in w 



rich an unstable dis 



short-lived clumps ( IPickett et al 



2003 



Mejfa et al 



2005 



£ frag ments into 



Durisen et al 



2008), no clumps 



ever form during the course of this simulation. Instead, the initial 2-, 3- and 4-armed 
perturbations grow briefly and then decay. Although the nonaxisymmetric structure reaches 
a modest nonlinear level, with amplitudes in 5p/p ~ 0.1 for m = 2, 3 and 4 at about 0.25 
ORP, no clumps form. Local cooling times t coo i for vertical columns through the disk in 
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units of the local disk orbital period P ro t vary from about 1 5 to 40 across t he most active 
regions of the disk (radial distances of 10 AU or less). The iGammid (120011 ) fragmentation 



criterion, as applied to a 7 = 5/3 EOS in a 3D simula tion, leads us to 



must be less than one for these regions to fragment (IRice et al. 



2003 . 



expec t that t coo y/Pj 



rot 



20051 ). Seeing no 



clump formation is consistent with our high values of t COQ \/P 1 



rot • 



We also in vestigated whether our disk is convectively unstable, because convection has 



been proposed (IBoss 



2004 



Mayer et al. 



20071 ) as an efficient cooling mechanism responsible 



for fragmentation. Our Schwarzschild criterion analysis shows that only small regions of 
the disk are susceptible to convective instability early in the evolution. By 1 ORP, these 
regions are confined to thin strips not associated with the dynamics of material near the 
equatorial plane. 



4. Discussion 

Our simulation results are clearly different from those highlighted in B07. All four 
simulations described in B07 have the same grid resolution (512 azimuthal cells) as ours. 
The B07 models H and TZ also use the same 50 K thermal bath BC. While other parameters 
are varied in B07 (for example, the condition of a free vertical temperature gradient or one 
forced to be monotonically decreasing), the variation in thermal treatment of disk models 
in B07 does not greatly affect the ultimate outcome of clump formation. 

As Boss notes in B07 and further explores, there are a number of possible reasons for 
the differences between his simulations and ours. Here, we have attempted to remove most 
of them. For example, because B07 did not implement an artificial viscosity, we turned 
off CHYMERA's artificial viscosity for the simulation reported here. We used the same 
disk, the same EOS, and the same initial perturbation as in B07. We also include envelope 
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irradiation to mimic the B07 thermal bath in a physically realistic manner. Although our 
code uses fewer terms in the boundary potential expansion than in B07, studies discuss ed 



in several of our paper s (IPickett et al. 



Boley fc Durisen 



2003 



Pickett fc Durisen 



2007 : 



Boley et al 



2007 



20081 ) show that our boundary potential expansion is adequate and that 



the 512 az imuthal cells used here are sufficient to detect fragmentation when it should occur 



^see, esp. 



Boley fc Durisen 



20081 ) . Moreover, the num erical heating in the interior disk 



regions that affected some of our previous simulations (jCai et al. 



(Cai et al. 


2006; 


Bolev et al. 


2006) 



has been eliminated by the increased effective vertical resolution of the current model, in 



which we typically resolve the vertical disk structure with at least 
and with many more near and outside 10 AU. Tests in iBoley et al. 



ce 



Is inside 10 AU 



( 120071 ) indicate that the 



disk must be vertically resolved by at least about 5 or 6 cells above the midplane to avoid 
numerical difficulties, including possible artificial heating. 

Given the similarity of the models and input physics and given the elimination of all 
other major numerical concerns, the radical disagreement in results between our simulation 
and those in B07 must be due to differences in treatment o f the photospheric BCs and 
the optically thin layers above the photosphere. iBossi (120041 ) recognizes that rapid cooling 



is n eeded for fragmenta tion and attributes this to thermal convection, which he and 



also 



Mayer et al. 



(120071 ) identify with upwellings at close to the sound speed associated 



with the spiral arms. The upwellings are actuall y hydraulic jumps due to disruption of 



vertical hydrostatic equilibrium in spiral shocks fBo lev fc Durisen 



produce rapid cooling when properly modeled (IRafikov 



2007 



2006) and should not 



Boley fc Durisen 



Boley et al. 



2006 



2007 



20081 ). In B07, however, the temperature excess over the thermal bath of 
any shock-heated material moving upward across the disk photosphere is effectively set to 
zero. To test how this affects cooling, we measured the flux of excess thermal energy in our 
post-shock regions at the time when we see maximum spiral wave amplitudes (top panel of 
Figure 1). We estimate that, in the densest shock structures, removal of this excess thermal 
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energy as the gas jumps above the photosphere would produce an effective local t coo \/P TO t 
of 2 to 3, compared with the true radiative t coo \/P Iot of more like 15 at the same places, a 
reduction in local cooling time by a factor of between 5 and 8! While t coo y/P TOt is not yet 
small enough to induce fragmentation, the B07 simulations do evolve denser structure and 
stronger shocks than our simulation beyond this time, and it is easy to imagine that the 
Gammie fragmentation criterion could then be satisfied. This estimate for the excess cooling 
introduced by the B07 temperature reset in the optically thin regions is highly suggestive 
that fast cooling in the B07 scheme is artificial. Unfortunately, when we tried to implement 
the B07 BCs in our own code, we encountered numerical instability in the radiative heating 
and cooling terms near the photosphere, and we could not test this hypothesis any further. 

So w hat then is t h e stat us of disk instability as a gas giant formation mechanism? 



Although 



Mayer et al. 



(120071 ) agree with Boss that disk instability may work under som e 



conditions inside 4 AU, simulations by oth er groups ( 



Forgan et al. 



20091 ) and analytic treatments (IRafikov 



Stamatc 



2005 



los k, Whitworth 



2008 



20071 ) support our arguments 



against fragmentation in the inner disks of young solar-type stars . On the other hand 



recent hydrodynamics simulat i ons (IStamatel 



2009 



Vorobyov fc B asu 



Rafikov 



2009 



2010 



Havfield et al 



Dodson-Robinson et al. 



os fc Whitworth 



2009 



Boley 



2009; 



20101 ) and other arguments ( 



Clarke 



Bol ev et al. 

T 



200C 



20091 ) now suggest that disk instability may very well 
operate in the outer regions of large (> 100 AU) massive disks during the early accretion 
phase to produce super- Jupiters, brown dwarfs, and low-mass stars. We think these 
fragmentation results in outer massive disks are quite plausible and worthy of further study. 
The disk instability mechanism that Boss has championed since 1997 may yet prove to be 
an important formation channel for gas giants, but in a different place than he originally 
proposed. 
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5. Summary 

We have conducted a high resolution simulation of a 4 to 20 AU protoplanetary disk 
of 0.091 M Q orbiting a 1 M Q star similar in most respects to simulations highlighted in 
B07. We have eliminated as many differences in methodology as possible, and we have 
answered the main arguments in B07 against the robustness of our previous results. The 
simulations in B07 show fragmentation into dense clumps; the simulation we present here 
shows no fragmentation. Because everything else is the same, the disagreement must be due 
to differences in the treatment of the optically thin regions and the photospheric boundary 
conditions. We identify at least one possible source of artificially fast cooling in the B07 
scheme and demonstrate that it can decrease cooling times by a large factor. This work 
supports the conclusion that gas giant planet formation will not occur via disk instability 
in the inner tens of AU of disks around solar-type stars. 
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Figure Captions 



Figure 1. Equatorial plane mass density contours for our disk simulation, in a format 
similar to B07. The contours span approximately 7 orders of magnitude in mass density, 
using factor of two spacing. The box that encloses each image spans 40 AU on each side, 
accommodating the initial disk radius of 20 AU. Shown are density images near the time 
of the highest amplitude nonaxisymmetry (top, 0.25 ORP), the time at the end of the 
equivalent runs in B07 (middle, 3.7 ORP), and the final image (bottom, 5.0 ORP). As in 
B07, we indicate the highest densities (> 10~ 10 gcm~ 3 ) with filled black contours. 
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